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Positional Disorder (Random Gaussian Phase Shifts) in the Fully Frustrated 
Josephson Junction Array (2D XY Model) 
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We consider the effect of positional disorder on a Josephson junction array with an applied magnetic 
field of / = 1/2 flux quantum per unit cell. This is equivalent to the problem of random Gaussian 
phase shifts in the fully frustrated 2D XY model. Using simple analytical arguments and numerical 
simulations, we present evidence that the ground state vortex lattice of the pure model becomes 
disordered, in the thermodynamic limit, by any finite amount of positional disorder. 
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The stability of vortex lattices to random disorder is a 
topic of considerable recent interest, motivated by studies 
of the high temperature superconductors. In two dimen- 
sions (2D), periodic arrays of Josephson junctions form 
a well controlled system for investigating similar issues 
of vortex fluctuations and disorder. Here we consider 
the effect of "positional" disorder on the vortex lattice of 
the fully frustrated Josephson array, with / = 1/2 flux 
quantum of applied magnetic field per unit cell. 

Positional disorder |l|-|j was first discussed with re- 
spect to the Kosterlitz-Thouless (KT) transition for the 
/ = model in zero magnetic field. Early arguments 
predicting a reentrant normal phase at low temperatures 
have been revised by recent works ^,|| which argue that 
there is a finite critical disorder strength a c ~ s/tt/8; 
for a < <j c an ordered state persists for < T < T c (a). 
For the pure / = 1/2 case on a square grid the 
ordered state has two broken symmetries: the U(l) sym- 
metry ("KT-like" order) associated with superconduct- 
ing phase coherence, and the Z(2) symmetry ( "Ising-like" 
order) associated with the "checkerboard" vortex lattice, 
in which a vortex sits on every other site. Previous works 
|0,|J have considered the effect of positional disorder on 
this / = 1/2 model; all have concluded that both Ising- 
like and KT-like order persist for at least small disorder 
strengths a. In this work, however, we present new argu- 
ments that suggest that, for / = 1/2, the critical disorder 
is <7 C = 0. 

The Hamiltonian for the Josephson array is given by 
the "frustrated" 2D XY model 101, 



H[0i] = U (°i - - A m) 

ifi 



(1) 



where i are the sites of a periodic square grid with basis 
vectors (jl — x, y, the sum is over all nearest neighbor 



(n.n.) bonds 



jl) and 6 1 — 9i+p. — A lfl is the gauge 



invariant phase difference across the bond, with = 
(</>o/2-7r) A • d£ the integral of the vector potential. 

Positional disorder arises from random geometric dis- 
tortions of the bonds of the grid, resulting in, — 



A 



(0) 



SAip,; A*fJ is the value in the absence of disorder, 



and 5Aifj_ is the random deviation. We take the SA^ to 



be independent Gaussian random variables with 
[SA ifl ] = 0, and [SA^SAjy] = a 2 8 l j5nv . 



(2) 



[. . .] denotes an average over the quenched disorder. The 
positionally disordered array is thus also referred to as 
the XY model with random Gaussian phase shifts. 

When U(4>) is the Villain function ||, the Hamiltonian 
(0) is equivalent to a dual "Coulomb gas" of interacting 
vortices HMO, 

H[th] = \ - / - 5fi)Gij{nj - / - 6fi) . (3) 

The sum is over all pairs of dual sites i, j, rii is the 
integer vorticity on site i, and the interaction Gij is the 
Green's function for the 2D discrete Laplacian operator, 
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-2Tr5ij, where Ay 
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ASij. For large separations, Gi 
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The fi = f 



Sfi are (l/27r) times the circulation of the 
/ is the average applied flux, while 
Sfi is the deviation due to the random SAi^, 



Aifj, around dual site i 



6f< - I ISA. 
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SA, 



SA, 
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Geometrically distorting a bond increases the flux 
through the cell on one side of the bond, while reduc- 
ing the flux through the cell on the opposite side by the 
same amount. The Sfi are thus anticorrelated among 
n.n. sites. Positional disorder is thus the same as ran- 
dom dipole pairs of quenched charges ±Sfi Q. From 
Eqs. (|) and M) we get, 



[Sfi\=0, and [SfiSfj] 



(5) 



The Hamiltonian (|3|) can be rewritten as interacting 
charges in a one body random potential S, 



n 



(6) 



where qi = m — f, and the random potential is Vi 
Ej Giitfi- For / = 1/2, 1i = ±1/2- From Eq. (§, 
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Wi]=0, and [V i Vj]=Y,Gik[5fkSfi]G l . 
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The Vi thus have logarithmic long range correlations. 

We now use an Imry-Ma |l2| type argument to esti- 
mate the stability of the doubly degenerate checkerboard 
ground state to the formation of a square domain of side 
L. The energy of such an excitation consists of a domain 
wall term, Ed, which is present for the pure case, and a 
pinning term, E p , due to the interaction with the random 
V t . E d {L) has the form 0, 



Ed ~ aL + c In L 



(8) 



The first term is the interfacial tension of the domain 
wall; the second term comes from net charge that builds 
up at the corners of the domain @. Calculating E d (L) 
numerically for a pure system, we find an excellent fit to 
Eq. (§), with a = 0.28, c = 0.15, and d = 0.058. 

By Eq. (Q), the average pinning energy of the domain 
V, [E p ] — 2^ igX ,<7i[Vi] = 0, but the variance is, 



charges qi , and the long range correlations of the random 
potential Vi. 

To check this prediction, we carry out Monte Carlo 
(MC) simulations of the Hamiltonian (0) with periodic 
boundary conditions on L x L square grids. Our MC 
procedure is as follows p5[ . One MC excitation attempt 
consists of the insertion of a neutral n = ±1 vortex pair 
on n.n. or next n.n. sites, which is accepted or rejected 
using the usual Metropolis algorithm. L 2 such attempts 
we call one MC pass. At each temperature we typically 
used 4000 MC passes to equilibrate the system, followed 
by 128, 000 MC passes to compute averages. Every 100 
passes we attempt a global excitation reversing the sign 
of all the charges, qi — > — q t . For each disorder realization 
we cooled down two distinct "replicas" , starting with dif- 
ferent random charge configurations and using different 
random number sequences. In only about 3% of the cases 
did the two replicas fail to give reasonable agreement. 

To test for Ising-like order we define an order parame- 
ter analogously to an Ising antiferromagnet, 
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(11) 
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where Eq = (tt/32)L 2 is the ground state energy of the 
checkerboard domain [ jl5| . The root mean square pinning 
energy is thus, 



[Ep^ rms — bL, 



b = 



2^2 



0.35cr 



(10) 



For domains whose energy is lowered by the interaction 
with Vi, the typical excitation energy is E = Ed — [E p ] rms . 
Eqs. d|) and ( |l0| ) imply that when b > a, i.e. when a > 
o~ c — 0.8, E(L) has a maximum at L = £ = (c/2^/2) /{a — 
a~ c ). Domains of size L > £ will lower their energy by 
increasing in size, and so disorder the system. Thus, one 
naively expects that when a < o~ c the system preserves its 
Ising-like order, but when a > a c the system is disordered 
into domains of typical size £. 

However the leading size dependencies of Eqs. (|8|) and 
(|l0|), Ed ~ [E p ] rms ~ L, are exactly the same as found 
in the 2D n.n. random field Ising model (RFIM). For the 
RFIM it is known 
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17] ] that 2D is the lower critical 
dimension, that the randomness causes domains walls at 
T = always to roughen and so acquire an effective 
negative line tension, and that the critical disorder is <7 C = 
0, i.e. any amount of disorder, no matter how weak, 
destroys the Ising-like order of the pure case. By analogy, 
we suggest that the positionally disordered / = 1/2 2D 
XY model similarly has a c = 0. Our conclusion, that 
\Ep\rms ~ L as in the 2D RFIM, follows from a subtle 
cancellation between the long range interactions between 
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FIG. 1. (a) [{M 2 }] vs. T, for a = 0.3, and sizes L = 10, 14, 
20 (solid lines are guides to the eye); Histogram of occurrences 
of (M 2 ) in 200 realizations of disorder, for a = 0.3 at low 
T = 0.02, for (b) L = 10, (c) L = 14 and (d) L = 20. 

We first consider a — 0.3, smaller than both the naive 
estimate of a c = 0.8, and the a c = ~ 0.63 of the 

/ = model. Fig. | plots [(M 2 )] vs. T, averaged over 
200 disorder realizations, for sizes L = 10, 14 and 20. All 
curves start to increase from zero near T ~ 0.13, which 
is T c (a = 0) of the pure model. However [(M 2 )] at low T 
decreases steadily with increasing L. The reason for this 
becomes clearer if we consider the histogram of values of 
(M 2 ) that occur as we sample the different realizations 
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of disorder. We show such histograms in Figs. [j]b-d, for 
the lowest temperature T = 0.02. As L increases, the 
statistical weight shifts from predominantly ordered sys- 
tems (M 2 = 1/4), to predominantly disordered systems 
(M 2 = 0). Assuming that this trend continues, we expect 
that as L -» oo, [(M 2 )} -> 0. 

To measure the "random field correlation length" £, we 
consider the vortex correlation function, 



1 



5 ( k ) = ^E elk ' (ri ~ rj) M 

i,3 



(12) 



For the pure case, S'(k) in the ordered phase has singular 
Bragg peaks at K = ±irx ± Try. If the vortex lattice 
is disordered, these peaks will broaden, and their finite 
width provides a measure of £. Writing k = K + <5k, and 
assuming a Lorentzian shape for the disorder averaged 
peak, [S'(k)] cx l/(5k 2 + £~ 2 ), we determine £ by fitting 
to this form for 5k = 0, and 5k = 2ir/L Q. In Fig. | we 
show £ vs. fj at our lowest T — 0.02, for several system 
sizes L. Only for our smallest value a — 0.25 does a finite 
size effect remain. In this case, however, £ decreases as L 
increases. This is in contrast to the increase of £ with L 
that one would expect if one were approaching a second 
order transition. This behavior is consistent with that 
seen in Figs. [j]b-d, where as L increases, a greater fraction 
of the disorder realizations result in disordered states. 
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FIG. 2. Correlation length f(cr) vs. a at T = 0.02 for 
various sizes L. Sizes L = 14, 20 are averaged over 200 real- 
izations of the randomness; sizes L = 30, 40 are averaged over 
50 realizations. Dashed, dotted and solid lines are fits to the 
scaling forms (i), (ii) and (iii) respectively. 

We next fit our results for £(<r) to several possible scal- 
ing expressions: (i) £ ~ e c ^ ct ~' Tc ^ , (ii) £ ~ e c/{<?-<?c) ^ 
and (iii) £ ~ \a — <J C \~ P - The first has been suggested 
by Binder §H§ for the 2D RFIM. While in Binder's ex- 
pression <t c = 0, here we leave it as an arbitrary pa- 
rameter to be determined from the fit. The second has 
been suggested for the positionally disordered / = 
model [0J|], in which er c > 0. The third is the famil- 
iar power law form. Using data for only the largest L 
for each cr, the results of these fits are shown in Fig. [|. 
The value of a c and the \ 2 01 the fit for each case is (i) 



cr c = 0.0046 ± 0.050, x 2 = 67; (ii) o c = 0.0134 ± 0.055, 
X 2 = 67; (iii) a c = 0.0013 ± 0.098, p = 2.86 ± 0.84, 
X 2 = 7.6. The power law (iii) gives a significantly better 
fit than (i) or (ii) , however all give a c = within the esti- 
mated error. Given the rather limited range of the data, 
the above fits should be treated with caution. However 
they do indicate that the data contains no suggestion of 
a diverging £ at a finite a. Coupled with our Imry-Ma 
argument, we thus find a consistent picture suggesting 
that cr c = for the / = 1/2 2D XY model. 
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FIG. 3. (a) Overlaps [Q se ifi], [Qseif2], and [<3 cross ] vs. T 
for a = 0.3 and L = 20; (b) Overlap susceptibility \Q vs - T 
for a — 0.3 and L — 10, 14, 20. Both are averaged over 200 
disorder realizations. 

Returning to the case a = 0.3, where Ising-like order 
has been lost, we now consider whether the system may 
still have a finite temperature "spin glass" transition to 
a disordered but frozen vortex state. To test for this we 
measure the self and cross overlaps JTsj ] , Qscif and Q C ross, 



(13) 



a and (3 index the two independent replicas. For r suf- 
ficiently large we expect Q se ifi = Q sc if2 = Qcross, if the 
system is well equilibrated. Averaging Eq. (|l^) over sev- 
eral values of r > 2000 to improve our statistics, we plot 
[Qselfl], [<9scif2] and [<2 cross ] vs. T in Fig. |a. We see 
that our system is fairly well equilibrated down to the 
lowest T we study. To test for a spin glass transition, we 
measure the overlap susceptibility, 



XQ=L 2 {{(Q 



2 ) 

cross/ 



[(Across) ] } 



(14) 



which we plot vs. T in Fig. |3jb for various system sizes. 
The peak in xq near T ~ 0.06 shows no noticeable in- 
crease as L increases, thus suggesting that there is no 
finite temperature spin glass transition. 

If the vortices are not frozen, but are free to diffuse, 
one expects that superconducting phase coherence is also 
destroyed. To explicitly test this we measure the helicity 
modulus. The Hamiltonian (||) can viewed as represent- 
ing the XY model with "fluctuating twist" boundary con- 
ditions . Using the method of Ref. EG] , we determine 
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the dependence of the total free energy F of the corre- 
sponding XY model, as a function of the twist (A x , A y ) 
which is applied in a "fixed twist" boundary condition. 
We then determine the (A^o, A y o) that minimizes F; the 
helicity modulus tensor is then the curvature of F at the 
minimizing twist, Y M „ = d 2 F/dA ll dA l/ . In Fig. [|a we 
plot Ti, the largest of the two eigenvalues of T M „, vs. T, 
for cr = 0.3 and sizes L = 10, 14, 20. At all T, Ti contin- 
ues to decrease as L increases, giving no suggestion of a 
finite temperature transition. In Figs, ^b-d we plot his- 
tograms of the minimizing twist Ao for the three sizes L. 
Note, in choosing our random phase shifts SA^, we im- 
pose the constraint JY SA^ — in order to remove one 
trivial source of A ^ 0. We see that the width of the 
distributions of Ao steadily increases with increasing L, 
suggesting [Q that the strength of the random disorder 
is renormalizing to greater values on larger length scales. 
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FIG. 4. (a) Ti vs. T, for a = 0.3, and sizes L = 10, 14, 20 
(solid lines are guides to the eye); Histogram of values of Ao 
found in 200 disorder realizations, for a = 0.3 at low T = 0.02, 
at (b) L = 10, (c) L = 14, and (d) L = 20. 

To conclude, our results suggest that Ising-like order 
is destroyed for any finite amount of positional disor- 
der. Further, we found in one specific case that when the 
Ising-like order vanished, no spin glass order or phase co- 
herence existed either. We speculate that this remains 
true as well for any finite disorder strength. Although 
<r c = 0, the finite £(cr) nevertheless can become extremely 
large for small values of a. When £ exceeds the size of 
the experimental or numerical sample, the system will 
indeed look ordered. We believe this explains previous 
numerical work on this problem which reported the per- 
sistence of Ising-like order at small cr. In the most recent 
of these works, Cataudella @] reports at cr ~ 0.113 a fi- 
nite T c to an Ising-like ordered state. The correlation 
length exponent that he finds is v ~ 1.7, clearly different 



from that of the pure model. Using our scaling form (iii) 
we can estimate that at this value of cr, £ ~ 120, much 
larger than Cataudella's largest system size of L = 36. 
His results may thus be reflecting a cross over region at 
L < £, rather than a true transition. 
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